Batch PRs #173, #175, #156, #158, #152 - #212
Merged
Merged
Conversation
`--soloUMIfiltering MultiGeneUMI_CR` kept every gene tied at the highest
read count. CellRanger's rule is the opposite on exactly that case: the
gene with the *strictly* highest count takes the UMI, and a tie means no
gene counts it.
STAR walks the genes keeping a running maximum and clears its winner
whenever it meets an equal count
(`SoloFeature_collapseUMIall.cpp:212-224`):
if (ig.second>maxu) { maxu=ig.second; maxg=ig.first; }
else if (ig.second==maxu) { maxg=-1; };
...
if ( maxg+1==0 ) continue; // not counted for any gene
One read per gene is the ordinary shape of a multi-gene UMI, and it is
always a tie, so the old rule made the flag inert in practice rather
than merely inaccurate. Measured on a 20 000-read 10x fixture (200 cells
from the real v3 whitelist, 400 genes, 720 UMIs deliberately shared
between two genes), against STAR 2.7.11b with the same flags:
identical entries STAR counts rustar counts
before 13 749 / 14 806 15 423 16 465
after 13 902 / 13 967 15 423 15 414
The flag removed nothing at all before; STAR removes 1 030 counts. The
gap goes from +1 042 to -9.
The outcome does not depend on the order the genes are visited — a
strict maximum always ends as the winner, a tie always ends with none —
so iterating a `HashMap` here stays deterministic.
`multi_gene_umi_cr_drops_a_tie_entirely` pins the case the old tests
missed: they only covered 3 reads against 1, where both rules agree.
Not yet implemented, and stated so rather than left to be discovered:
STAR applies a second condition, that the winning gene must also hold
the top count among *uncorrected* UMIs (`umiGeneMapCount0`, same file,
lines 226-232). That needs the pre-correction counts, which this code
does not keep. The 65 entries still differing out of 13 967 are the
place to look for its effect.
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
STAR corrects UMIs within each gene *before* deciding which gene owns a UMI, and applies two conditions, not one (`SoloFeature_collapseUMIall.cpp:134-148` and `:203-235`): 1. one gene must hold a strictly higher read count than every other, on the **corrected** UMI map — that is #173, already landed; 2. and that winner must not be beaten in the **uncorrected** map at the same key. The second condition exists because correction moves reads between UMIs: a gene can win only because correction folded a neighbouring UMI onto it, and STAR rejects that win rather than counting it. Reproducing it needs the order STAR uses. The generic path here filters multi-gene UMIs first and corrects afterwards, which cannot express either condition: by the time correction happens the ownership decision is already made. `MultiGeneUMI_CR` therefore takes its own path, which is also what STAR does — the flag is only valid with `--soloUMIdedup 1MM_CR`, so there is no combination this bypasses. `cellranger_1mm_map` exposes the correction mapping that `cellranger_1mm` already computed and threw away. Measured against **CellRanger 10.0.0** on the 20 000-read fixture from #172, with #165 and #173 also applied: identical entries CellRanger rustar #165 + #173 13 651 / 13 709 15 111 15 091 plus this change 13 676 / 13 709 15 111 15 116 Entries CellRanger has and we do not go from 29 to 7, and the count gap from -20 to +5, which is 0.03%. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
…_distribution STARsolo's `EmptyDrops_CR` rescue draws from `std::mt19937`, converts to doubles with `std::generate_canonical<double, 53>`, and picks categories with `std::discrete_distribution`. Two of those three are implementation-defined in the parts that matter: the standard fixes mt19937's output but not how `generate_canonical` consumes it, and says nothing about how `discrete_distribution` maps a uniform onto categories. So porting "the algorithm" is not enough — it has to be libc++'s algorithm, because that is what STAR is built against and where its numbers come from. libc++ accumulates two 32-bit draws in *ascending* significance and divides by 2^64; a most-significant-first accumulation, or one draw scaled to 53 bits, both give perfectly good uniforms and neither reproduces STAR. Every expected value in the tests came out of a C++ program compiled against the real libc++ and run, not from reading its source. `tests/libcxx_oracle.cpp` is that program, kept so the values can be regenerated rather than trusted. `generate_canonical` is compared as bit patterns, since a difference in the last place changes which category a sample lands in. Not yet wired into the EmptyDrops path. `solo::count` samples with a `SplitMix64` stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree on an arbitrary number of draws. It is true of whatever was checked and unknown elsewhere. Replacing it moves cell calls, so it belongs in its own change with the solo differential run against it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Two approximations in the CellRanger cell-calling path are replaced by what CellRanger and STAR actually compute. Both move cell calls, which is the point: the previous numbers were plausible rather than right. The ambient profile is now smoothed with Simple Good-Turing (Gadsby & Sampson, via Elworthy's implementation, which is what STAR vendors). The ambient counts come from a small sample of empty droplets, so a gene seen twice there is not twice as likely as one seen once, and a gene seen zero times is not impossible — it is one the sample was too small to show. SGT fits the frequency spectrum and reserves mass for the unseen from the singleton rate, then smooths the rest along a log-log line. What was here before had the right shape and the wrong numbers: it reserved mass the same way but distributed the remainder in proportion to raw counts, with no smoothing at all. The Monte-Carlo null is now drawn with libc++'s `std::mt19937` and `std::discrete_distribution`, seeded `19760110 * (isim + 1)` per simulation, as STAR seeds it. The previous sampler was a SplitMix64 stream under a comment calling it "WeightedIndex-equivalent; empirically byte-identical EmptyDrops cell calls" — a claim that cannot hold in general, since two unrelated generators cannot agree over an arbitrary number of draws. The libc++ types were ported and checked against real libc++ in the previous commit on this branch; this wires them in. One generator per simulation, no shared state, so the walks still run in any order on any number of threads and give the same p-values. D17 comes with it: STAR leaves `PZero` uninitialised when the spectrum has fewer than five distinct frequencies and `analyse()` bails, so it reads whatever the stack held. Here it is zero from construction, which is what "no basis for reserving unseen mass" means. Recorded in docs-old/dev/divergences.md.
Section 1.2, in the What STAR does / What rustar-aligner does / Why / Impact / Source format CONTRIBUTING.md asks for, replacing the docs-old file the earlier version of this work carried.
Quantifies transcripts rather than genes, from where each read's 3' end sits
relative to each transcript's. In a 3'-biased assay that distance is what
separates isoforms: a read 200 bases from the end of one and 4000 from the end
of another is evidence for the first. The distribution of those distances is
estimated from the run's own histogram, smoothed and cut where the 3' peak
decays into the body, and used as the likelihood in an EM over UMIs.
Concordance needed no new code. `align_to_transcripts` already refuses to
project an alignment that leaves the transcript, touches an intron, or crosses a
junction the transcript does not have — which is exactly STAR's `Concordant`
(`Transcriptome_classifyAlign.cpp`). A projection that survives is concordant;
one that does not, is not. The projection also puts the 5' end at coordinate
zero for both strands, so the distance to the 3' end is one expression rather
than two.
Two behaviours worth stating because they are not the obvious ones:
Output is per cluster, not per cell, and `--soloClusterCBfile` is required.
A single cell does not have enough UMIs to resolve isoforms, so the EM would be
fitting noise. Asking for the feature without a clustering is refused rather
than run.
A UMI seen on several reads contributes the *intersection* of their transcript
sets. Those reads came from one molecule, so a transcript missing from any of
them cannot be its source. Taking the union would let a single stray read
resurrect an isoform every other read excluded.
Two of STAR's quirks are reproduced rather than corrected, because the cut point
and every weight depend on them: the running-average divisor is `min(2N+1,
i + N)` rather than the number of elements actually summed, and the transcript
length factor is taken from the cumulative distribution at `trLen - 1`
(`SoloFeature_quantTranscript.cpp`).
Numbers are formatted the way C++'s default stream prints them — six
significant digits, fixed inside `[1e-4, 1e6)` and scientific outside — since
the normalised distribution runs down to ~1e-4 where Rust's `{}` and C++'s
default disagree on both notation and digit count.
Records are accumulated under a mutex, so there are no partials to merge. It was dead from the moment it was written; CONTRIBUTING.md rules out shipping it.
parametersDefault puts both Transcript3p and --soloClusterCBfile between "#####UnderDevelopment_begin : not supported - do not use" and "#####UnderDevelopment_end", and STAR --help prints that banner around them. The module said none of this. It matters for how the port is read: it follows STAR's code, so it inherits the unfinished parts of that code, and a differential against STAR compares two implementations of something STAR does not support. A reviewer should be told that before deciding to take it.
…GeneUMI `--soloUMIfiltering MultiGeneUMI_All` resolved to the same variant as `MultiGeneUMI`, which is neither what STAR does nor what the option is documented to do. Of the three available behaviours it was the only one nobody had asked for. In STAR the option is a no-op: it is parsed and stored, but its consumption site tests only the `MultiGeneUMI` flag, so selecting it leaves the filter entirely off. Documented, it removes a UMI seen in more than one gene from *all* of them, rather than from the losers only. `UmiFiltering::MultiGeneUmiAll` now exists and does the documented thing: a UMI appearing in several genes is evidence of a collision or of chimeric amplification, so it is discarded outright rather than attributed to whichever gene happened to read deepest. Single-gene UMIs are untouched, which the test checks across every mode. Raised upstream as #144 before changing it, since "be faithful to STAR" and "do what the flag says" genuinely point in opposite directions here. Also adds `docs-old/dev/divergences.md`, recording this and the homopolymer-UMI rule, so deliberate differences are written down rather than rediscovered as surprises in a differential run. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Batch merge of the solo/UMI-counting PRs from @BenjaminDEMAILLE, plus the fixes found while validating them. Each of his commits is preserved with original authorship; the batch is split by commit so it reads PR by PR, with fixes as separate commits.
Closes #173, #175, #156, #158, #152.
What's in it
MultiGeneUMI_CRkept every gene tied at the highest read count; CellRanger gives a tied UMI to no gene. One read per gene is the ordinary shape of a multi-gene UMI, so the flag removed nothing in practice.EmptyDrops_CRon CellRanger's actual statistics: Simple Good-Turing ambient profile, and a bit-exact libc++mt19937/discrete_distributionfor the Monte-Carlo null.--soloFeatures Transcript3p, with--soloClusterCBfile.MultiGeneUMI_Alldoes what it documents instead of aliasing toMultiGeneUMI.Changes made while merging
A missing validation (new). STAR refuses
--soloUMIfiltering MultiGeneUMI_CRunless the dedup is exactly1MM_CR(ParametersSolo.cpp:463-468); we accepted it silently and counted with whatever map the chosen dedup produced. This is item 3 of #172, and it belongs here because the batch is what makes the flag do anything — before #173 it was inert, so the missing guard could not mislead anyone. We now emit STAR's message verbatim.Four DIVERGENCE.md corrections. The entries arriving with these PRs needed work:
MultiGeneUMI_All's entry said STAR's flag is unimplemented. It is implemented, and matches the documentation exactly (SoloFeature_collapseUMIall.cpp:79-88zeroes every gene for a UMI seen in more than one). What fails is the site that acts on those counts::116gates onMultiGeneUMI, and the two flags are set in mutually exclusive branches (ParametersSolo.cpp:457-462), so selecting_Allzeroes the counts and never reads them. That makes this a one-line wiring bug rather than a design decision — we are diverging from STAR's behaviour but not from its intent. Recorded with the line numbers.1.2. Filed under section 1 and renumbered.3.2citedmerging_partials_preserves_order_and_sums, removed by solo: --soloCellReadStats CB writes CellReads.stats #157's own cleanup in Batch PRs #166, #170, #171, #157, #159 #200. Dropped the citation; the entry still holds.Not taken: the
keep only this PR's entrycommits. Each stacked PR carries one, stripping the other PRs' CHANGELOG or DIVERGENCE entries — correct on his branch, destructive in a batch. #158's would have deleted #156's CHANGELOG entry. #152's strips #165's divergence entry, which we do want gone since #165 is not here, so that one was applied by hand instead.Validation
Default paths byte-identical to
main. Records compared with headers excluded, since the@PG CL:line necessarily differs by output prefix:mainAgainst STAR 2.7.11b, unchanged from the recorded baselines:
The new guard refuses the same command STAR refuses, with the same message.
Gate: 617 tests,
cargo clippy --all-targetsclean,cargo fmt --checkclean.What we could not validate
#175 and #156 are unmeasurable on the data we have. Both are byte-identical before and after on the 10x mouse chr19 set:
MultiGeneUMI_All,MultiGeneUMI_CRand the default all produce 379 entries / 433 UMIs, because no UMI in this dataset is shared between genes.So the claim for these three is faithful port, no regression — not "validated". Each was instead checked line by line against STAR's C++, which is recorded in the commits.
#172 item 4 is the fix, so it would be great to get a reproducible method for it from @BenjaminDEMAILLE (over guessing)
One assumption this invalidates
The single EmptyDrops cell we disagree with STAR on has long been explained as Monte-Carlo noise. #156 makes the sampler bit-exact with STAR's and that cell still disagrees, so the explanation is probably wrong. I'll look into this in a separate PR.